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USING  THE  DWOPER  ROUTING  MODEL 
TO  SIMULATE  RIVER  FLOWS  WITH  ICE 

Steven  F.  Daly  and  George  D.  Ashton 

INTRODUCTION 

The  National  Weather  Service  is  responsible  for  forecasting  water 
flows  and  water  levels  of  the  major  rivers  of  the  United  States.  Part  of 
the  forecast  methodology  involves  calculation  of  the  expected  water  levels 
and  flows  based  on  existing  flows.  To  perform  the  calculations  a  variety 
of  numerical  models  are  used,  all  of  which  are  based  on  the  principles  of 
open  channel  hydraulics.  The  most  sophisticated  of  the  currently  used 
models  is  DWOPER  (Dynamic  Wave  Operational  Forecast  Program)  (Fread  1976). 
However,  none  of  the  models  include  the  effects  of  ice  on  the  hydraulics  of 
the  river  flows  and  consequently  on  the  water  levels  and  discharges  to  be 
expected  when  ice  is  present  in  a  river  The  effort  reported  here  is  a 
first  attempt  to  include  the  effects  of  ice  on  the  predicted  water  levels 
and  flows  in  a  model  that  includes  unsteady  (i.e.  time-varying)  effects. 

DWOPER  is  an  implicit  dynaml  .-uting  model  that  is  being  used  by  the 
National  Weather  Service  to  forecast  floods  and  other  flow  conditions  on 
major  wacerways.  The  model,  which  predicts  flows  and  stages,  is  based  on 
the  one-dimensional  differential  equations  of  unsteady  flow  known  as  the 
Saint-Venant  equations.  These  equations  consist  of  an  equation  of  continu¬ 
ity,  which  conserves  the  mass  of  the  flow,  and  an  equation  of  momentum, 
which  conserves  the  flow  momentum.  Unlike  earlier,  simpler  forecast 
models,  DWOPER  uses  the  complete  unsteady  flow  equations  which  include  the 
effects  of  frictional  resistance,  flow  accelerations  and  surface  slope. 
DWOPER  also  provides  additional  hydraulic  information  about  flow  along  the 
waterways  such  as  depth,  flow  cross-sectional  area  and  top  width,  hydraulic 
radius,  velocity,  water  surface  slope,  and  energy  slope. 


DWOPER  does  not,  however,  include  the  effects  of  ice  on  the  hydraulics 
of  the  flow.  River  ice  can  have  a  large  effect  on  the  stages  and  flows  of 
a  river.  In  fact,  every  ice  process  that  occurs  on  a  river  produces  a 


hydraulic  transient  of  some  magnitude.  When  an  ice  cover  is  present  on  a 
river,  the  frictional  resistance  of  the  cover  to  the  flow  must  be  taken 
into  account.  A  certain  portion  of  the  channel  cross  section  will  also  be 
blocked  by  the  ice  cover  and  be  unavailable  for  flow.  In  the  case  of 
severe  ice  jams  this  portion  may  be  quite  large.  The  hydraulic  radius, 
which  is  defined  as  the  cross-sectional  area  of  the  flow  divided  by  the 
wetted  perimeter,  must  be  modified  to  account  for  the  increase  in  the 
wetted  perimeter  caused  by  the  ice  cover. 

Also,  the  more  "dynamic"  ice  effects  should  be  taken  into  account. 
These  include  ice  cover  initiation,  ice  cover  growth  and  evolution,  ice 
jamming,  ice  cover  breakup  and  ice  cover  movement.  The  hydraulic  tran¬ 
sients  associated  with  these  events  can  be  very  large.  The  conditions 
leading  to  the  initial  bridging  of  a  river  by  floating  ice  are  not  well 
known,  unless  there  is  an  obvious  surface  obstruction  such  as  a  floating 
ice  boom  or  a  dam.  Factors  that  influence  bridging  include  ice  discharge 
rate,  fragment  size,  open  water  width,  ice  concentration,  lateral  ice 
growth,  air  temperature  and  the  ice  transport  capabilities  of  a  given  river 
reach. 

On  the  other  hand,  the  mechanics  of  the  subsequent  ice  cover 
accumulation  process  have  been  studied  to  a  much  greater  extent.  Once  an 
ice  jam  or  bridging  has  begun,  the  ice  cover  accumulates  through  the 
arrival  of  ice  fragments  or  blocks  at  the  upstream  end  of  the  ice  cover. 

The  ability  of  the  arriving  blocks  to  resist  entrainment  in  the  water 
flowing  under  the  existing  ice  can  be  estimated  by  knowledge  of  the  thick¬ 
ness  of  the  blocks  and  the  velocity  of  the  water.  At  higher  velocities  the 
blocks  will  underturn  and  be  swept  under  the  ice  cover;  they  will  deposit 
and  thicken  the  cover  somewhere  downstream.  The  "mechanical  stability"  of 
the  cover  to  withstand  the  forces  applied  it  by  the  flowing  water,  gravity 
and  wind  can  also  be  estimated.  If  the  cover  is  not  strong  enough,  it  will 
collapse  and  thicken  until  it  has  enough  internal  strength  to  resist  the 
applied  forces. 

Once  an  initial  cover  has  been  formed,  further  thickening  results  from 
loss  of  heat  to  the  atmosphere.  This  thermal  ice  growth  can  be  estimated 


from  air  temperature  and  wind  speeds  by  use  of  relatively  simple  formula¬ 
tions. 

Ice  cover  breakup  is  a  process  that  is  not  understood.  Thermal  and 
mechanical  effects  both  play  large,  if  only  qualitatively  understood, 
roles  in  it. 

To  summarize,  river  ice  and  the  hydraulic  features  of  .iver  flow  form 
an  imperfectly  understood  feedback  loop.  The  hydraulic  characteristics  of 
the  flow  affect  the  formation  and  evolution  of  the  ice  cover,  which  in  turn 
affect  the  hydraulic  characteristics  of  the  flow.  Ice  effects  on  river 
flow  can  be  quite  dramatic,  and  serious  efforts  are  being  made  to 
understand  and  quantify  these  effects  and  their  consequences.  It  is  the 
purpose  of  this  paper  to  present  a  method  of  including  the  effects  of  river 
ice  in  DWOPER,  to  take  some  of  the  initial  steps  needed  to  include  those 
effects,  to  show  examples  of  the  associated  hydraulic  transients,  and 
finally  to  suggest  additional  steps  that  can  be  taken  and  indicate  where 
further  research  is  required. 

CHANGES  NEEDED  TO  ADAPT  DWOPER  TO  ICE  CONDITIONS 
Hydraulic  parameters 

Two  basic  equations,  the  equation  for  conservation  of  the  mass  of  the 
flow  and  the  equation  for  conservation  of  the  momentum  of  the  flow, 
compose  the  one-dimensional  differential  equations  of  unsteady  flow,  known 
as  the  Saint-Venant  equations,  which  are  the  theoretical  bases  of  DWOPER. 

To  fully  account  for  the  presence  of  an  ice  cover  and  the  transport  of  the 
ice  in  the  channel,  both  equations  must  be  modified.  The  conservation  of 
mass  equation  must  be  modified  to  not  only  conserve  the  mass  of  the  water 
flow,  but  also  the  mass  of  the  ice  flow.  The  conservation  of  momentum 
equation  must  conserve  the  momentum  and  forces  acting  on  the  ice  cover  when 
it  is  moving,  as  well  as  the  momentum  and  forces  acting  on  the  flow. 
However,  at  this  stage  of  model  development,  we  will  assumed  the  ice  cover 
to  be  either  stationary  or  nonexistent.  That  is,  the  influence  of  the 
motion  of  the  ice  will  not  be  included  in  either  equation.  Therefore,  we 
will  examine  the  conservation  of  momentum  equation  and  only  the  modifica¬ 
tions  required  to  account  for  the  presence  of  a  stationary  ice  cover. 

The  law  of  conservation  of  momentum  is  used  to  derive  the  second 
Saint-Venant  equation,  i.e.,  the  equation  of  motion  or  dynamic 
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PROFILE 


Figure  1.  Definition  sketch  of 
control  volume  for  derivation 
of  momentum  equation  (from  Fread 
1976). 

equilibrium.  The  conservation  of  motion  is  given  by  Newton's  second  law  of 
motion,  which  may  be  stated  as:  The  sum  of  the  forces  acting  on  the 
surface  of  the  control  volume  plus  the  net  rate  of  momentum  entering  the 
control  volume  equals  the  time  rate  of  accumulation  of  momentum  within  the 
control  volume. 

The  forces  acting  on  the  surface  of  the  control  volume  of  length  dx 
shown  in  Figure  1  include  1)  the  gravity  force  due  to  the  weight  of  the 
fluid  Fg,  2)  the  force  due  to  the  frictional  resistance  along  the  channel 
bottom  and  sides  Ff,  3)  the  force  due  to  the  shear  stress  produced  by 
wind  movement  at  the  free  surface  of  the  control  volume  Fw,  and  4)  the 
unbalanced  pressure  force  Fpj  -  Fpr;  $  is  the  bottom  slope  angle  and 
FpW  Iz  the  pressure  force  on  the  sides  due  to  enlargement  or  constric¬ 
tion. 

Of  the  four  forces  accounted  for  in  the  above  description  only  the 
force  due  to  the  frictional  resistance  along  the  channel  is  directly 
affected  by  a  stationary  ice  cover  floating  in  the  control  volume.  A 
stationary  cover  would  also  interfere  with  the  wind  force.  However,  the 
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wind  force  is  usually  either  negligible  or  unknown  and  will  not  be  included 
in  this  discussion. 

The  frictional  resistance  is  manifested  by  a  shear  stress  t  along  the 
bottom  and  sides  of  the  control  volume.  An  empirical  equation  for  open 
channel  resistance,  the  Manning  equation,  is  used  to  express  the  frictional 
resistance.  The  Manning  equation  (English  units)  is 

v  .  hJHX  e2/3  s1/2  (1> 

n.  r 


where  V  is  the  average  velocity  over  the  flow  cross  section,  R  is  the 
hydraulic  radius,  Sf  is  the  friction  slope,  and  n^  is  the  Manning 
roughness  coefficient  of  the  channel  bottom  with  dimensions  of  ft^/^. 
The  conservation  form  of  the  momentum  equation  is 
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where 

Q  *  flow  rate 

t  *  time 

3  *  momentum  correction  coefficient 

x  =  distance  along  the  channel 

A  *  area  of  the  flow  cross  section 

g  =  gravity 

y  ■  depth  of  flow 

SQ  =  bottom  slope 

q  *  lateral  inflow  rate  with  a  velocity  component  Vx  in  the  x 
direction 

tw  =  wind  shear  stress  at  the  water  surface 
B  =  width  of  flow 

p  =  water  density. 

The  friction  slope  Sf  can  be  evaluated  by  the  Manning  equation. 
Substituting  Q  =  VA  into  eq  1,  we  find 
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The  hydraulic  radius,  which  is  defined  as  A/P  where  P  is  the  wetted  perim¬ 
eter,  can  be  adequately  approximated  by  the  hydraulic  depth  H  =  A/B  for 
most  natural  river  channels  for  which  B  >  lOy.  Therefore,  we  can  write 


,  n2  Q2 
f  2.21  A2  H4/3 


However,  we  must  modify  A,  R  and  n  in  the  above  expression  to  account  for 
the  presence  of  the  ice  cover. 

The  cross-sectional  area  of  the  flow,  A,  must  be  reduced  to  account 
for  the  reduction  of  the  cross-sectional  area  by  the  ice  cover.  If  we 
assume  a  uniform  ice  cover  of  thickness  n,  the  flow  area  A^  under  the  ice 
is 


A  =  A - 1  nB  (5) 

i  P 

and  A  is  determined  as  the  cross-sectional  area  corresponding  to  the  water 
level  with  ice  present,  is  the  density  of  ice  and  B  is  the  width  of 
the  cross  section. 

The  hydraulic  radius  of  a  channel  is  defined  as  the  cross-sectional 
area  of  the  channel  divided  by  wetted  perimeter.  The  presence  of  an  ice 
cover  will  increase  the  wetted  perimeter  by  B,  the  width  of  the  cross 
section.  If  we  again  assume  that  the  depth  is  small  when  compared  to  the 
width,  the  hydraulic  radius  with  an  ice  cover  Rj  is 


The  Manning  roughness  coefficient  is  an  empirical  measurement  of 
channel  roughness.  It  must  be  modified,  relative  to  open  water  conditions, 
to  account  for  the  roughness  of  the  underside  of  the  ice  cover.  It  should 
be  emphasized  that  the  roughness  of  the  ice  cover  can  vary  both  in  space 
and  through  time.  It  can  consist  of  both  small  scale  (surface)  drag  and 
large  scale  drag  when  the  underside  is  a  jumbled  pile  of  angular  blocks  and 
slabs.  Many  methods  of  determining  the  composite  roughness,  nc,  of  an 
ice  covered  channel  have  been  proposed.  The  following  formula  is  easy  to 
use  and  produces  reasonable  results  as  shown  by  Carey  (1967)  and  Uzuner 
(1975): 
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(7) 
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where  is  the  roughness  coefficient  of  the  underside  of  the  ice  cover 
and  nj,  is  the  roughness  coefficient  of  the  channel  bottom. 

Equation  2,  the  momentum  equation,  can  now  be  rewritten  for  the  case 
of  an  ice-covered  river,  with  the  wind  stress  term  dropped,  as 
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and  eq  4  can  be  rewritten  as 
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Ice  parameters 

Initiation  of  cover 

As  was  noted  in  the  introduction,  unless  there  is  an  obvious  surface 
obstruction  such  as  a  floating  ice  boom  or  a  dam,  the  conditions  leading  to 
the  initial  bridging  of  a  river  by  floating  fragments  are  not  well  known. 
Factors  that  influence  bridging  include  ice  discharge  rate,  fragment  size, 
open  water  width,  ice  concentration,  lateral  ice  growth,  air  temperature, 
and  the  ice  transport  capabilities  of  a  given  river  reach.  Experience  has 
shown  that  on  given  rivers  the  ice  cover  will  initiate  at  about  the  same 
places  every  winter  season.  These  sites  are  generally  those  places  where 
the  water  surface  slope  changes  from  steeper  to  milder.  The  reasons  for 
this  are  not  entirely  clear  at  this  time;  however,  to  successfully  model  a 
given  phenomenon,  clear  and  unambiguous  criteria  must  be  provided.  Three 
approaches  are  proposed  to  achieve  this:  an  arbitrary  initiation,  use  of 
the  Pariset  and  Hausser  stability  criteria,  and  a  criterion  based  on 
surface  concentration  of  floating  ice. 

Arbitrary  initiation.  A  site  where  an  ice  cover  will  begin  can  be 
chosen  arbitrarily.  The  choice  can  be  based  on  experience,  judgment  or 
investigative  need.  It  is  a  simple  programming  matter  to  select  a  location 
in  the  waterway  where  the  ice  cover  will  be  initiated. 
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Pariset  and  Hausser  stability  criteria.  In  a  now  classic  paper 
Pariset  et  al.  (1966)  developed  criteria  to  determine  the  thickness  of  an 
ice  cover  formed  by  accumulation  of  floating  ice  fragments.  Their  analysis 
presumed  that  bridging  had  already  taken  place  and  examined  the  thickness 
of  an  ice  cover  that  progresses  upstream  by  accumulation.  The  essence  of 
the  analysis  is  a  balance  between  the  water  forces  on  the  floating  ice 
cover  and  the  resisting  forces  between  the  ice  cover  and  the  banks,  with 
the  thickness  determined  by  the  internal  strength  of  the  accumulation  that 
is  necessary  to  transmit  the  forces  to  the  banks.  While  their  criteria 
were  developed  for  the  case  of  ice  cover  progression  rather  than  initia¬ 
tion,  it  may  be  argued  that  moving  ice  of  the  same  thickness  as  the  ice 
that  they  were  describing  will  bridge  across  a  channel  and  stop  since  that 
ice  has  the  potential  of  resisting  the  forces.  The  criteria  are  examined 
in  more  detail  later.  For  the  moment,  we  simply  note  that  the  criteria 
relate  the  water  depth,  water  discharge  Q,  channel  width  B  and  the 
hydraulic  roughness  C  to  a  function  of  the  ratio  of  the  ice  cover  thickness 
n  to  depth  of  flow  H  in  the  form 

x  =  BC^iF  =  f(ri/H)  *  (10) 

Pariset  and  Hausser  used  the  Chezy  coefficient  C  for  the  roughness  param¬ 
eter  associated  with  the  ice  undersurface  but  it  may  easily  be  related  to 
the  Manning  coefficient  (n  =  R^^/C). 

The  function  f(n/H)  is  shown  in  Figure  2.  If  its  value  falls  within 
or  on  the  boundaries  of  the  stable  region,  a  cover  is  considered  to  be 
stable  and  initiation  is  presumed  to  have  taken  place.  Experience  gained 
during  the  present  study  showed  this  to  be  a  very  difficult  criterion  to 
satisfy;  moving  ice  doesn't  often  fall  within  the  stable  region. 

Concentration  criteria.  As  ice  first  forms  and  moves  down  a  river,  it 
floes  together  at  the  surface  forming  ice  pans  and  hence  is  confined  to  the 
upper  layer  of  the  flow.  If  it  remains  only  one  layer  thick,  the  ice  must 
satisfy  a  continuity  of  surface  flux  of  ice  from  one  point  to  another  rela¬ 
tion.  Expressing  the  surface  concentration  of  ice  as  c,  the  surface  veloc¬ 
ity  as  V  and  the  width  as  B,  then  at  points  1,2  we  find 

cx  Vl  Bj  *  c2  V2  B2  .  (11) 
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Figure  2.  Stability  diagram 
(after  Pariset  et  al.  1966). 


With  knowledge  of  the  variation  in  V  and  B  along  a  reach,  the  concentration 
of  ice  of  any  point  may  be  determined  if  the  concentration  is  known  at  some 
upstream  point.  If  this  concentration  is  greater  than  the  critical  concen¬ 
tration  at  which  bridging  takes  place,  we  then  may  presume  that  the  ice 
cover  has  been  initiated. 

It  is  difficult,  however,  to  determine  the  initial  ice  concentrations 
and  the  value  of  the  critical  concentration.  Limited  work  has  been  done  on 
the  latter  problem.  Frankenstein  and  Assur  (1972)  speculated  that  broken 
ice  will  clog  and  come  to  a  standstill  when  the  surface  concentration  is 
about  0.7,  but  they  offered  no  data  to  support  this  value.  More  recently 
Ackermann  et  al.  (1979)  analyzed  the  ice  transport  capacity  of  a  river. 

They  included  the  interaction  of  ice  floes  and  were  able  to  calculate  the 
limiting  capacity  of  a  given  channel  to  transport  ice.  In  all  cases 
examined  the  maximum  ice  discharge  corresponded  to  concentrations  of  0.6  < 
c  <  0.7  (Ackermann  et  al.,  in  press).  They  also  found  by  numerical  analy¬ 
sis  that  the  maximum  Ice  discharge  depended  on  the  characteristic  diameter 
of  floes  (larger  floes  clog  more  easily)  and  the  slope  (the  milder  the 
slope  the  less  the  required  ice  discharge).  The  most  troublesome  parameter 
for  our  purposes  here  is  the  floe  diameter  and  this  brings  us  full  circle, 
i.e.,  even  if  criteria  based  on  ice  concentration  (or  floe  diameter)  are 
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available,  it  is  currently  beyond  the  capability  of  river  ice  mechanics  to 
determine  these  parameters,  other  than  by  direct  observation. 


Ice  cover  progression 


Once  an  ice  cover  has  initially  bridged  a  channel  the  cover  may  grow 
upstream.  The  stability  of  ice  blocks  arriving  at  the  upstream  end  of  the 
cover  can  be  assessed  using  a  simple  equilibrium  analysis.  If  the  velocity 
of  the  approaching  block  is  above  some  critical  velocity  Vc,  the  block 
will  underturn  and  be  sw.’pt  under  the  cover;  if  below  the  critical  veloc¬ 
ity,  the  block  will  come  to  rest  against  the  upstream  end  of  the  cover  and 
the  cover  will  progress  upstream.  Ashton  (1974)  found  the  critical 
velocity  to  be 


V 

c 


2(i  -r) 

nT  2  1/2 

[5  -  3(1  -  ] 


(12) 


where  pj  and  p  are  the  densities  of  the  ice  and  water  respectively,  and 
ni  is  the  thickness  of  the  ice  block.  If  the  approach  velocity  is  above 
the  critical  velocity  for  underturning,  we  can  estimate  the  resultant 
thickness  n  of  the  jam  by  solving  the  expression  determined  by  Pariset  and 
Hausser  (1961), 


P,  1/2 

V  "  [2g(l  -  (l-§).  03) 


Rates  of  ice  cover  growth  can  then  be  simply  determined  by  balancing  the 
incoming  ice  flux  with  the  quantities  required  to  increase  the  cover  at 
this  thickness. 

Thickening 

The  ice  cover  may  initially  thicken  further  by  mechanical  collapse  in 
response  to  the  increasing  forces  of  shear  stress  on  the  underside  of  the 
Initial  accumulation.  Once  this  process  is  completed,  the  ice  cover 
further  thickens  by  thermal  growth. 

As  the  ice  cover  progresses  upstream  by  accumulation  at  the  thickness 
n ,  which  is  found  by  application  of  eq  13  using  the  local  variables  of 
velocity  and  depth,  it  may  not  have  enough  internal  strength  to  transmit 


Figure  3.  Definition  sketch 
for  analysis  of  ice  accumulation 
stability. 
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the  shear  stresses  acting  on  its  underside  to  the  banks.  This  is  particu¬ 
larly  true  for  so-called  "wide"  rivers.  If  the  ice  cover  is  not  strong 
enough  it  will  thicken  further  by  a  process  often  termed  "shoving"  until  it 
has  the  competence  to  transmit  the  forces  to  the  banks.  Pariset  and 
Hausser  (1961)  analyzed  this  situation  as  follows. 

A  differential  length  dL  (Fig.  3)  of  an  ice  cover  of  thickness  n  and 
width  B  experiences  an  overall  force  balance  of 

Bdf  +  2  (x  n  +  pf)dL  =  (x  +  x  )BdL  (14) 

c  ^  ^ 

where  tc  is  the  cohesion  of  the  ice  to  the  banks,  p  is  a  type  of  ice 
friction  coefficient  on  the  banks,  xw  is  the  water  shear  stress  on  the 
underside  of  the  cover,  and  xg  is  the  gravity  force  per  unit  area  acting 
on  the  cover  along  the  slope.  Equation  14  may  be  integrated  to  yield 


f 


B(x  +  x  ) 

r  w  g 

[  2p 


x  n 


exp(-  -^)] 


(15) 


where  the  boundary  condition  that  f  *  0  at  L  *  0  has  been  used. 

A  distinction  can  be  made  in  eq  15  between  wide  and  narrow  rivers.  If 
the  first  bracketed  term  is  negative  the  river  is  characterized  as  narrow 
because  the  ice  cover  can  carry  the  shearing  forces  of  the  water  and  the 
gravitational  forces  without  further  thickening. 

In  eq  14  and  15,  xw  can  be  calculated  from  open  channel  hydraulics, 
as  can  xg.  Pariset  et  al.  (1966)  found  from  field  measurements  that  p 
had  a  value  of  about  1.28  and  the  product  xcn  varied  from  75  to  92  lb 
ft-1  for  the  St.  Lawrence  River  (but  it  varied  for  other  rivers). 
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If  the  river  is  wide  the  cover  will  thicken  by  shoving  until  the 
internal  resistance  is  sufficient  to  carry  the  forces  to  the  banks. 

Once  an  ice  cover  has  been  initiated  on  a  river  and  is  not  moving, 
it  thickens  further  by  thermal  growth.  Incremental  thickening  can  be  quite 
well  estimated  by  stepped  integration  of 


dn  _  1  r 

dt 


T  ) 
m 


ia 


(16) 


where 

n  *  ice  thickness 
t  *  time 

p±  **  density  of  ice 
X  =  heat  of  fusion 

ns  =  thickness  of  any  snow  layer  on  the  surface 
k^  and  kg  =  thermal  conductivities  of  ice  and  snow 
Ta  =  air  temperature 
Tm  =  32° F,  the  melting  point  of  ice 
hia  «  heat  transfer  coefficient  applied  to  the  difference  between  the 
top  surface  temperature  and  the  air  temperature. 

Equation  16  is  based  on  a  quasi-steady  heat  conduction  analysis  (Ashton 
1980)  that  is  more  than  adequate,  considering  the  uncertainty  in  pragmatic¬ 
ally  available  values  of  the  input  parameters.  The  value  of  h^a  ideally 
is  calculated  using  detailed  energy  budget  analyses;  as  a  practical  matter 
hia  may  be  estimated  and  is  of  the  order  of  3  to  4  Btu  ft-2  hr-1  °F-1. 

HYDRAULIC  TRANSIENTS  ASSOCIATED  WITH  ICE  PROCESSES 
Introduction 

Hydraulic  transients  are  associated  with  all  ice  processes  that  occur 
on  a  river.  The  magnitude  of  the  associated  transients  can  be  very  large 
or  very  small,  depending  on  a  variety  of  conditions.  In  this  report  we 
investigate  the  potential  magnitude  of  the  hydraulic  transients  by  examin¬ 
ing  two  processes:  the  freezeup  and  the  breakup  of  a  uniform  ice  cover.  An 
idealized  prismatic  channel,  with  a  variable  slope  and  hydraulic  roughness, 
is  used  to  examine  the  hydraulic  transients  associated  with  these  two 
processes.  One  of  the  difficulties  of  studying  the  unsteady  effects  of  ice 
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development  is  the  almost  complete  lack  of  reliable  field  data.  Therefore, 
our  first  purpose  is  to  determine  what  is  important  and  what  is  not. 

An  important  concept  in  the  study  of  transients  caused  by  an  ice  cover 
is  the  concept  of  normal  or  uniform  depth.  The  depth  of  a  uniform  flow  is 
the  normal  depth.  A  given  river  discharge  will  be  at  normal  depth  when  the 
slope  of  the  water  surface,  Sf ,  is  equal  to  the  slope  of  the  channel 
bottom;  therefore 


Sf  -  s0 

where  SQ  is  the  slope  of  the  channel  bottom, 
uniform  depth,  we  know  from  eq  3 

^2/3  ‘/2 

n  o 


(17) 


For  open  water  flow  at 


(18) 


and  that  for  wide  channels 
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where  HN  is  the  normal  depth.  Rearranging,  we  find 
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(19) 


(20) 


Similarly,  for  an  ice  covered  river,  assuming  that  the  hydraulic  roughness 
coefficient  of  the  ice  undersurface  is  the  same  as  the  bottom  and  that  we 
can  ignore  temporarily  the  portion  of  the  channel  blocked  by  ice,  we  find 


or 


q  -  Hr* 

“nice  ■  u32hn  • 


(21) 


(22) 


Thus ,  the  presence  of  an  ice  cover  increases  the  normal  depth  by  approxi¬ 
mately  32%.  If  the  area  of  the  channel  blocked  by  the  ice  is  included 
(i.e.  the  depth  is  called  the  free  water  surface),  it  has  been  found  that 
an  accurate  description  for  the  normal  depth  is 


HNICE 
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where  p^n/p  represents  the  displacement  effect  of  the  ice  cover. 

While  it  is  seldom  that  a  given  discharge  is  actually  flowing  in  a 
channel  at  its  normal  depth,  it  is  the  change  of  normal  depth  as  a  result 
of  an  ice  cover  forming  or  breaking  up  and  the  actual  water  depth  relative 
to  the  resultant  normal  depth  that  determines  the  potential  magnitude  of 
the  associated  hydraulic  transient.  The  ultimate  magnitude  of  the  hydraul¬ 
ic  transient,  however,  will  be  determined  by  the  downstream  constraints, 
the  slope  of  the  river  channel,  the  length  of  the  channel,  and  the  duration 
of  the  ice  processes.  It  is  hard  to  disentangle  the  relative  importance  of 
these  factors;  however,  the  channel  slope  is  probably  the  most  important 
factor.  The  steeper  a  channel,  the  shorter  will  be  the  distance  over  which 
the  transition  to  a  new  equilibrium  depth  occurs  in  response  to  a  change  in 
the  downstream  conditions  (an  example  is  shown  in  Fig.  4).  Hence,  for  a 
given  downstream  constraint,  the  slope  of  the  channel  largely  determines 
the  extent  of  influence  of  an  ice  formation  or  ice  breakup  event.  The 
larger  transients  are  expected  on  the  steeper  channels. 

Breakup 

Figure  5  illustrates  the  transient  response  of  three  prismatic  rec¬ 
tangular  channels  with  varying  channel  slopes  to  breakup  of  a  uniform  ice 
cover.  The  downstream  elevation  in  each  case  is  held  constant  at  the  ice 
covered  normal  depth,  and  the  upstream  discharge  held  constant.  The  ice 


Figure  4.  Example  calculation 
of  distance  upstream  at  which 
99%  of  normal  depth  with  ice 
is  attained  with  downstream 
constraint  fixed  at  open  water 
normal  depth 


Figure  5.  Transient  response  to  ice 
breakup  for  prismatic  channels  of  three 
different  slopes. 


Figure  6.  Effect  of  length  of  chan 
nel  on  transient  response  to  break¬ 
up  (S  =  0.001). 


cover  was  assumed  to  instantaneously  cease  being  a  hydraulic  influence. 

(In  the  numerical  simulation  this  was  done  by  instantaneously  changing  the 
hydraulic  radius  to  the  depth  of  flow  beneath  the  ice  cover,  reverting  from 
the  composite  roughness  to  the  more  usual  open  channel  roughness,  and 
adjusting  the  cross-sectional  area.)  Only  for  the  steeper  channels  does 
this  approximation  approach  reality.  The  magnitude  of  the  transient 
response  does  indeed  increase  with  channel  slope.  The  milder  channels  also 
tend  to  show  more  reflections  of  the  transient,  and  the  oscillations  of  the 
discharge  can  be  noted  for  these  cases.  The  effect  of  varying  the  length 
of  the  channel  on  the  transient  response  can  be  seen  in  Figure  6.  The 
larger  transient  is  due  to  the  cumulatively  longer  release  of  storage  water 
caused  by  the  change  in  resistance.  In  general,  the  longer  and  steeper  the 
channel  is,  the  greater  the  transient  response  to  ice  breakup. 

Freezeup 

We  can  reverse  the  above  process  to  study  the  transient  response  to 
the  freezeup  of  a  river  channel.  Figure  7  illustrates  the  transient 
response  of  the  same  three  prismatic  rectangular  channels  to  the  ice  cover 
freezeup.  The  downstream  elevation  in  each  case  is  held  constant  at  the 
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Figure  7.  Transient  response  to  instantaneous 
freezeup  of  rectangular  prismatic  channels  of  three 
different  slopes  (Q  »  300,000  ft3  s-1,  constant 
downstream  elevation  at  open  water  normal  depth). 

open  water  normal  depth,  and  the  upstream  discharge  held  constant.  The 
hydraulic  influence  of  the  ice  cover  was  assumed  to  appear  instantaneously 
over  the  entire  length  of  the  channel.  We  know  that  river  freezeup  is 
generally  a  gradual  process  that  happens  over  a  length  of  time.  Thus, 
these  examples  are  extreme  cases  to  examine  the  potential  magnitude  of  the 
hydraulic  transients  that  may  take  place.  Immediately  apparent  Is  the 
numerical  instability  shown  for  the  steepest  river.  This  will  be  more 
fully  discussed  in  the  next  section.  Again,  the  magnitude  of  the  transient 
response  increases  with  channel  slope.  Many  different  processes  can 
combine  to  cause  a  river  to  freezeup.  The  river  velocity  it  known  to  be  a 
critical  parameter,  and  the  lateral  growth  rates,  ice  inflows  from  up¬ 
stream,  and  meteorological  conditions  all  contribute  to  the  overall 
process.  The  situation  is  greatly  simplified  in  this  simulation,  but 
again,  it  allows  examination  of  potential  magnitudes. 

MODIFICATIONS  TO  DWOPER  PROGRAM 

The  modifications  to  DWOPER  are  very  straightforward.  Basically,  the 
terms  of  the  conservation  of  momentum  equation  are  modified  as  described 
earlier.  The  hydraulic  influence  of  the  ice  cover  is  "felt"  by  the  channel 
through  the  modification  of  three  parameters:  the  composite  hydraulic 
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Figure  8.  Effects  of  weighting  factor  on  stability 
for  a  case  of  instantaneous  freezeup  (S  =  0.001). 


roughness  of  channel,  the  cross-sectional  area  and  the  hydraulic  radius. 
Whenever  the  presence  of  an  ice  cover  was  indicated  by  an  ice  thickness 
greater  than  zero,  these  parameters  were  changed  to  reflect  that  presence. 

Several  runs  were  made  to  test  the  effect  of  varying  the  weighting 
factor,  9.  This  parameter  determines  the  positioning  of  spatial  deriva¬ 
tives  and  nonderivative  terms  between  adjacent  time  lines.  Generally,  a 
value  for  0  of  0.55  to  0.6  is  used,  and  this  was  the  value  used  in  the  pre¬ 
vious  examples.  However,  an  instability  problem  was  noted  for  the  steepest 
channel  during  freezeup.  No  instability  problems  were  noted  for  the  case 
of  ice  breakup.  Problems  of  instability  of  the  numerical  solutions  are  a 
result  of  the  nonlinearity  of  the  various  terms  of  the  governing  equa¬ 
tions.  Frictional  resistance  is  a  very  nonlinear  term  that  assumes  more 
importance  for  the  steeper  channels.  The  presence  of  an  ice  cover 
increases  frictional  resistance.  Therefore,  the  freezeup  on  a  steeper 
channel  would  be  a  most  likely  process  to  cause  instability  problems  to 
appear.  Figure  8  shows  the  results  of  several  runs  that  were  made  with 
varying  values  of  0.  We  can  see  that  the  instabilities  did  not  entirely 
disappear  until  0  was  set  equal  to  a  value  of  1.00.  However,  the  front  of 
the  wave  was  slightly  lagged  and  the  magnitude  damped.  This  would  suggest 
that  a  variable  0  that  can  change  in  anticipation  of  or  in  response  to  the 
changing  ice  conditions  may  be  appropriate. 
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SUMMARY 

This  report  is  a  step  toward  an  eventual  goal  of  accurately  modeling 
the  processes  of  river  ice  and  the  associated  transient  responses.  We  have 
shown  that  by  modifying  the  equation  of  the  conservation  of  the  momentum  of 
the  flow  by  modifying  three  parameters  —  the  hydraulic  roughness,  cross- 
sectional  area  and  the  hydraulic  radius  —  we  can  indeed  produce  hydraulic 
transients  that  qualitatively,  at  least,  resemble  the  real  world.  We  are  a 
long  way,  however,  from  starting  our  simulation  in  the  fall,  and  predicting 
the  ice  processes  and  resultant  transients  through  to  the  spring  breakup 
period.  The  following  points  summarize  our  findings. 

1.  The  transient  response  of  channels  to  the  process  of  ice  freezeup 
and  ice  breakup  can  be  demonstrated  by  modifying  the  appropriate  parameters 
in  the  equation  of  the  conservation  of  momentum  of  the  flow.  In  general, 
the  steeper  the  channel,  the  greater  the  transient  response. 

2.  A  more  fully  developed  and  robust  simulation  of  the  transient 
response  could  probably  be  achieved  by  including  the  conservation  of  ice 
discharge  and  the  conservation  of  ice  momentum  and  forces.  A  variable  0, 
the  weighting  factor,  that  could  be  varied  in  response  to  the  changing  ice 
conditions  would  be  appropriate. 

3.  To  accurately  simulate  the  ice  cover  growth,  the  physics  and 
mechanics  of  lateral  ice  growth,  longitudinal  ice  floe  growth,  ice  cover 
initiation  by  bridging  or  arching,  and  ice  cover  breakup  must  be  more  fully 
understood. 

4.  The  analytical  investigation  of  the  transient  response  of  channels 
and  rivers  to  ice  processes  must  be  broadened  and  deepened.  Very  little 
work  has  been  done  in  this  area. 

5.  The  collection  of  field  data  must  be  improved  because  very  little 
field  data  are  available.  While  the  situation  is  improving,  the  lack  of 
methods  for  remote  sensing  of  ice  thicknesses  is  severely  hampering  the 
field  effort. 

6.  The  effect  of  ice  jams,  as  opposed  to  uniform  ice  cover,  should  be 
investigated.  Hydraulically,  a  jam  can  be  defined  as  any  ice  cover  that 
produces  a  steep  gradient  relative  to  the  channel  slope  across  the  length 
of  the  cover.  The  breakup  of  jams  approaches  the  dam  break  problem  as  a 
wor at  case.  The  rapid  formation  of  jams  produces  a  response  that  closely 
resembles  that  produced  by  partial  closure  of  a  sluice  gate  as  a  worst 
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case;  positive  waves  propagating  in  the  upstream  direction  and  negative 

waves  propagating  downstream  from  a  sudden  jam  site  have  been  measured  in 

the  field. 
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